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ABSTRACT 

Because of mass assignment onto grid points in the measurement of the 
power spectrum using the Fast Fourier Transform (FFT), the raw power spec- 
trum (\o~f{k)\ 2 ) estimated with FFT is not the same as the true power spec- 
trum P(k). In this paper, we derive the formula which relates (\S^(k)\ 2 ) to 
P{k). For a sample of N discrete objects, the formula reads: {\5^{k)\ 2 ) = 
J2 n [\W(k+2k N n)\ 2 P(k+2k N n) + l/N\W(k+2k N n)\ 2 }, where W(k) is the Fourier 
transform of the mass assignment function W(r), k^ is the Nyquist wavenumber, 
and n is an integer vector. The formula is different from that in some of previous 
works where the summation over n is neglected. For the NGP, CIC and TSC 
assignment functions, we show that the shot noise term Y2 n 1 / N\W (k + 2kN'n)\ 2 ] 
can be expressed by simple analytical functions. To reconstruct P(k) from the 
alias sum J2 n |W / (k + 2/cjvn)| 2 P(k + 2/cjvn), we propose an iterative method. We 
test the method by applying it to an N-body simulation sample, and show that 
the method can successfully recover P(k). The discussion is further generalized 
to samples with observational selection effects. 

Subject headings: methods: data analysis - methods: statistical - galaxies: clus- 
tering - large-scale structure of Universe 

1. INTRODUCTION 

The power spectrum P(k) of the spatial cosmic density distribution is an important 
quantity in galaxy formation theories. P(k) on large scales is a direct measure of the pri- 
mordial density fluctuation, and P(k) on small scales carries information of later non-linear 
evolution, therefore measuring P(k) can serve to distinguish between different theoretical 
models. 
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The power spectrum P(k) as a clustering measure has already been applied by many 
authors to observational samples of galaxies and of clusters of galaxies, including the Cf A and 
Perseus-Pisces redshift surveys (Baumgart & Fry 1991), the radio galaxy survey (Peacock 
& Nicholson 1991), the IRAS QDOT survey (Kaiser 1991), the galaxy distribution in the 
nearby superclusters (Gramann & Einasto 1992), the SSRS redshift sample (Park, Gott, & 
da Costa 1992), the CfA and SSRS extensions (Vogeley et al. 1992; da Costa et al. 1994), 
the 2Jy IRAS survey (Jing & Valdarnini 1993), the 1.2Jy IRAS survey (Fisher et al. 1993) 
and the redshift samples of Abell clusters (Jing & Valdarnini 1993; Peacock & West 1992). 
The power spectrum is also widely measured for cosmological N-body simulations, since it 
can easily characterize the linear and non-linear evolutions of the density perturbation (e.g. 
Davis et al. 1985). 

All of these workers except Fisher et al., have used the Fast Fourier Transformation 
(FFT) technique to make the Fourier Transforms (FT). In fact, one can use the direct 
summation [Eq. (5) below] to measure the power spectrum for a sample of a few thousand 
objects (Fisher et al. 1993). However it appears to be more convenient for most workers to 
use FFT packages (now available on many computers) to analyze the power spectrum. This 
is probably the reason why most of the previous statistical studies of the power spectrum 
have used FFT. For N-body simulations, one has to use FFT to obtain the power spectrum, 
since normally there are more than a million particles. When using FFT, one needs to 
collect density values p(r 9 ) on a grid from a density field p(r) (or a particle distribution), 
which is usually called 'mass assignment'. The mass assignment is equivalent to convolving 
the density field by a given assignment function W(r) and sampling the convolved density 
field on a finite number of grid points. The FT of p(r 9 ) generally is not equal to the FT of 
p(r), and the power spectrum estimated directly from the FT of p(r g ) is a biased one. The 
smoothing effect has already been considered in many previous works (e.g. Baumgart & Fry 
1991; Jing & Valdarnini 1993; Scoccimarro et al. 1998), but the sampling effect has not. 

In this paper, we derive the formulae which express these effects of the mass assignment 
on the estimated power spectrum. A procedure is proposed to correct for these effects, in 
order to recover the true power spectrum. The procedure is tested and shown to be very 
successful with a simulation sample. 

2. Formulae 

Let us first recall the definition of the power spectrum P(k). Let p(r) be the cosmic 
density field and p the mean density. The density field can be expressed with a dimensionless 
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fleld 5(r) (which is usually called the density contrast): 

Based on the cosmological principle, one can imagine that p(r) is periodic in some large 
rectangular volume V^. The FT of 5(r) is then defined as: 

<f(k) = ±r I 5(r)e- k rfr, (2) 

and the power spectrum P(k) is simply related to 5 (It) by 

p (k) = (I S(k) | 2 ) , (3) 
where (• • •) means the ensemble average. 

In practical measurement of P(k) either for an extragalactic catalogue or for simulation 
data, the continuous density field p(r) is sampled by only a finite number N of objects. In 
these cases, one needs to deal with the 'discreteness' effect arising from the Poisson shot 
noise. To show this, let's consider an ideal case where no selection effect has been introduced 
into the sample (In fact, some kind of selection effects must exist in extragalactic catalogues. 
We will discuss it in Section 3). The number density distribution of objects can be expressed 
as n(r) = J2j$ D ( r ~ r j)> where Vj is the coordinate of object j and 5 D (r) is the Dirac-5 
function. In analogy with the continuous case, the FT of n(r) is defined as: 

^(k) = ±[ n(r)e^dv - , (4) 

where n is the global mean number density, the superscript d represents the discrete case 
of p(r), and 5 K is the Kronecker delta. Following Peebles (1980, sections 36-41), we divide 
the volume into infinitesimal elements {dVj} with rii objects inside dVi. Then the above 
equation can be written as: 

^(^^E^-^. ( 5 ) 

i 

where N is nV M , the number of objects in V^. Since dVi is taken so small that rij is either 
or 1, we have rii = nf = nf = • • •, (rij) = ndVi and (niUj)^ = n 2 dVidVj[l + (5(ri)5(Tj))]. 
We can find the ensemble average of 5 d (k 1 )5 d *(k 2 ): 

(5^)5'*^)) = ^^(^,)e^- kl -^ k2 -^ 1)0 ^ 2)0 
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= (5(k 1 )5*(k 2 )) + ^ uk2 . (6) 
The last equation assumes ki ^ or k 2 ^ 0. The true power spectrum is then: 

p (k) = (I S(k) | 2 ) =(|^(k)| 2 )-^. (7) 

So the discreteness (or shot noise) effect is to introduce an additional term 1/N to the power 
spectrum (| 5 d (k) | 2 ). This fact is already well known to cosmologists. We present the above 
derivation because this method is useful in the following derivations. 

In principle one can use the direct summation of Eq. (5) to measure the power spectrum 
for a sample of discrete objects. However, as described in Section 1, most of the previous 
statistical studies of the power spectrum have used FFT. Moreover it would be impossible 
to use direct summation to measure the power spectrum for an N-body simulation. The 
quantity computed by the FFT is 

Sf ( k ) = ^J2 nf ^ etr9 ' k -<o, (8) 
g 

where the superscript / denotes quantities in FFT. n^(r g ) is the convolved density value on 
the g-th grid point r g = gH (g is an integer vector; H is the grid spacing): 

n f (r g ) = J n(r)W(r-r g )dr, (9) 

where W(r) is the mass assignment function. 

Following Hockney and Eastwood (1981), Eq.(8) can be expressed in a more compact 
way by using the so-called 'Sampling Function'. The sampling function II (r) is defined 
as a sum of the Dirac-5 functions spaced at unit length in all three spatial directions, i.e. 
n(r) = Z]{ n } 3° ( r — n )> where n is an integer vector. Defining 

n' f (r) = n(-|) J n( ri )H/(n - r)dn (10) 



and constructing 
one can easily prove 



^(k) = i- /n^(r)e- k rfr-^ , (11) 



N 

5 ,f (k) = 5 f (k) . (12) 
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Therefore one can express <^(k) as [cf. Eqs. (10-12)]: 

*'(k) = 1 f U^^mWin - r)e^dr - 8« . (13) 

The ensemble average of <5^(ki)<5^*(k 2 ) then reads: 
(^(k0^(k 2 )) = ± f n(|)n(|)[^(^)W(r il )W(r, 2 ) + ^(n,)W(r a )W(r i2 )" 

xe iri kl - ir2 k2 d ri dr 2 - 1 J n(^) ^(n l )^(r, 1 )e iri - kl 5£ i0 rfr 1 

/ n (|)E^)^(^) e " tr2 - k2 «,o^ 2 + «,oC )0 , (14) 

where = — r^-. Using 

n(k) = ±r f n(^)e ir k rfr = J2 <W> (15) 
where = k/H is the Nyquist wavenumber, one can find: 

(^(k^*^)) = £ [|W(k' 1 )| 2 P(k' 1 )^ lik , a + ^\W(k\)\X^ 2 ] , (16) 



ni,n2 



where k • = kj + 2k N ia. i and W(k) is the FT of W(r). For k x = k 2 = k we obtain our desired 
result: 

(|^(k)| 2 ) = \W(k + 2k N n)\ 2 P(k + 2k N n) + 1 £ |W(k + 2A^n)| 2 , (17) 

n n 

where the summation is over all 3D inetger vectors n. The meaning of the above equation is 
very clear. The density convolution introduces the factor VU 2 (k) both to the power spectrum 
and to the shot noise (1/iV). The finite sampling of the convolved density field results in 
the the 'alias' sums (i.e., the sums over n). The alias effect is well-known in Fourier theory, 
but has not been taken seriously in the power spectrum analysis of large-scale clustering in 
observational cosmology. Both effects of the convolution and the alias are significant near 
the Nyquist wavenumber kjy, see Figs. 1 & 2 in the next section. 



3. A PROCEDURE TO RECOVER P(K) 



In the practical measurement of P(k) using FFT, one should first choose a mass as- 
signment function. The NGP (p — 1), CIC (p = 2) and TSC (p = 3) assignment functions 
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are the most popular function for this purpose. For these schemes, we have (Hockney & 
Eastwood 1981): 



sin(^) sin(^) sinf^ 



/ 7rfci \ / wk 2 \ I irk 3 \ 
\2k N >\2k N >\2k N > 



WCk) = — ^ l kN . 2kN ' ' , (18) 

where ki is the i-th component of k. 



Once one has selected the assignment function, the shot-noise effect [the second term 
on the rhs of Eq.(17)] is easy to correct. For the NGP, CIC, and TSC assignments, the shot 
noise term can be expressed by: 

n 

where Ci(k) are simple analytical functions: 

( 1, NGP; 
C 1 (k)=\ n,[l-lsin 2 (^)], CIC; (20 ) 

lnai-sm 2 (^) + ^sin 4 (^)], TSC. 

Furthermore can easily find out that the Ci(k) of the CIC and TSC schemes are approxi- 
mately isotropic for k < fcjy, i.e. 

J [1- |sin 2 (^)], CIC; 
C M"{[l-^(£) + ±^(g-)], TSC. (21) 

We have tested Equations (19-21) by calculating (| 5^(k) | 2 ) for 10 random simulation 
samples, each of which consists of N = 10 5 points randomly distributed in a unit cube. In 
this case, -D 2 (k) = (| 5* (k) | 2 ). In the upper panel of Figure 1, we show the average values 
and la errors of (D 2 (k)N/Ci(k))d estimated from the 10 samples, where (and below) (■ • -) d 
means an average over all directions of k. In this calculation, we have used the NGP, CIC and 
TSC assignment functions, and have used Eq. (20) for Ci(k). The result for each assignment 
is shown by one symbol in the figure. From Eq. (19), one expects (D 2 (k)N/Ci(k)) d = 1 
for all the three assignment functions. Clearly the simulation results agree very well with 
Eq. (20). In the lower panel, we show the results in a slightly different way, i.e. the averages 
and la errors of (D 2 (k)N)d- The solid line is (D 2 (\t)N)d = 1 for the NGP assignment. The 
dotted and dashed lines are the approximate expressions of Eq. (21) for the CIC and TSC 
assignments. Again we find a very good agreement between Eq. (21) and the results from 
the random sample. This means that in most applications, one can use Eq. (21) to correct 
the shot noise in the FFT measurement of the power spectrum. 
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After correcting the shot noise, our central problem becomes how to extract P(k) from 
the first term of Eq. (17). Let us consider the correction factor C 2 (k) which is defined as 

(En ^ 2 ( k + 2k N n)P(k + 2k N n)\ 
C(k) = \ m U . ( 22 ) 

Since W(k) is a decreasing function, and P(k) is also a decreasing function on the scales 
k > k N , we expect that the alias contribution of large | n | to C 2 (k) is small. In particular, 
for k <C kN, any alias contribution is small and we have C 2 (k) tv W 2 (k) ~ 1 indpendent 
of P(k). For k ~ k^, the alias contribution to C2{k) becomes important, most of which is 
owing to the alias of |k + 2/cjvn| ~ k. Therefore the dependence of C^k) on P(k) is only the 
shape of P(k) at k ~ fcjv, i.e. the local slope c*n of P(k) at k ~ fcjv, a at = [lnP(fc)/ ln/c]^^. 

Since the local slope is unknown a priori in practical measurement of P{k), we 
propose an iterative method to get the correction factor C^/c). Suppose that we have 
measured the power spectrum P r (k): 

Pr(k) = (<|^(k)| 2 >-/J 2 (k)) d 

= (j2 w2 ( k + 2k Nn)P(k + 2k N n)^ (23) 

n 

for k < fcjv- The local slope cto of P r {k) at /c ~ A;^ is calculated by a power-law fitting to 
P r (/c) at 0.5/cjv < k < fcjv- Assuming a power-law form A; a ° for P(/c) of Eq. (20), we calculate 
C2(k,ao) and get Po(/c) = P r (k)/C2(k,ao). Using the local slope ao at 0.5/cjv < k < kN 
of the Po(k) just obtained, we calculate C 2 (k,a ) and Po(/c) = P r (k)/C 2 (k, a ) again. This 
calculation is repeated until Po(k) (or a ) converges to some defined accuracy. The converged 
Po(k) is our desired P(k). 

Success of the above iteration procedure is shown by Figure 2, where we present our 
measurement of P(k) for a simulation sample of Jing et al. (1995). The sample is five 
realizations of P 3 M simulation of a low-density flat universe with f2 = 0.2, A = 0.8 and 
h — 1. The simulation box size is 128/i -1 Mpc and the number of simulation particles is 
N = 64 3 . Since the details of the simulation are unimportant for our result here, we will 
not discuss them anymore. The upper panel of the figure shows the raw power spectra 
(|<5^(k)| 2 ), which are calculated with the NGP, CIC or TSC mass assignment and with 256 3 
or 64 3 grid points. As expected, the raw power spectrum depends on the scheme of the 
mass assignment: higher order mass assignment gives a smaller (|<5^(k)| 2 ) near the Nyquist 
frequency. At k ~ k^, the Nyquist wavenumber of 64 3 grid points, the raw power spectra 
calculated with 256 3 grid points are expected be very close to the true power spectrum P(k) 
because neither effect of the shot noise [P^ 2 (A;) P(k)], the convolution [W(fc) 1] or 
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the alias [C 2 (k) 1] is important. The differences at k ~ k^f between the power spectra 
calculated with 64 3 grid points and those with 256 3 grid points, show the importance of the 
effects discussed in this paper. The lower panel of the figure plots the power spectra P{k) 
which are corrected following the procedure prescribed previously in this section. In this 
example, we use Eq. (21) to correct the shot noise, and we require that the slope a® converge 
to the accuracy of |Aao| < 0.02 in the iterative method of deconvolving the alias summation. 
For this accuracy, only fewer than 5 iterations are needed in each calculation. The six power 
spectra P(k), measured with different mass assignments and with different numbers of grid 
points, agree so well that their curves overlay each other in the figure. The biggest difference 
between the 64 3 and 256 3 P(k) is at k = k^f, which is less than 4 percent. The result is very 
encouraging, and it tells us that using the correction procedure prescribed in this paper, one 
can obtain the true P{k) for k < k^ in the FFT measurement, independent of which mass 
assignment is used. 

In the above derivations and discussions, we have assumed, for simplicity, that the 
sample is uniformly (within Poisson fluctuation) constructed in a cubic volume. However, 
the above method is also valid for a sample with selection effects (e.g., an extragalactic 
catalogue). To show this, let us introduce a selection function S(r) which is defined as the 
observable rate of the sample at position r. If the underlying density distribution is n(r), 
the density distribution n s (r) of the sample is then: 

n s (r) = S(r)n(r) . (24) 

In this case, we define the following transformation for FFT: 

^( k ) = ^E^^)-^M e ^' k ' ( 25 ) 

9 

where n{(r g ) is the convolved density of n s (r) at r g [cf. Eq. (9)], and n is the mean underlying 
number density. Following the derivation of Eq. (16), we can easily find: 

(tfCkOtf'fka)) = ^[^(k' 1 )H/(-k' 2 )^5(k' 1 + k'3)P(k'3)5(-k' 2 -k' 3 ) 

ni,n 2 k' 3 

+ -U(k' 1 - k' 2 )W(k' 1 )W(-k' 2 )} , (26) 



N 



where S(k) is defined as: 



S(k) = JsWr J S(r)e "' k ' ir - (27) 

Because S(k) peaks at k 0, the coupling of the selection and the power spectrum [the first 
term of Eq. (26)] is important only at small k. How to treat this coupling is a nontrivial 
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task and beyond the scope of this paper (but see e.g. Peacock & Nicholson 1991; Jing & 
Valdarnini 1993). However for larger k where the effects discussed here become significant, 
we have: 

(^(k)! 2 ) » ^( k 3) Yl W ^ k + 2 ^n)P(k + 2^11) + ^ W ^ k + 2 ^ n ) ( 28 ) 

k 3 n n 

The difference between Eq. (26) and Eq. (19) is only the factor S 2 (ks). Therefore our 
procedure for correcting the effects of the mass assignment is still valid for observational 
samples with selection effects. Equation (28) is derived on the assumption that S(k) is 
compact in k space. This assumption is valid for many redshift surveys, e.g., the CfA and 
IRAS redshift surveys. But is not valid for surveys of irregular boundaries, e.g., pencil 
beam surveys. In addition, current redshift surveys of galaxies already cover a volume of 
(1000 /i _1 Mpc) 3 . Because one is also interested in the clustering information down to scales 
of 0.1 ft, _1 Mpc, a (10, 000) 3 gridpoint FFT is required to explore the clustering on all scales, 
which is still difficult to be realized on modern supercomputers. Other estimators, e.g., the 
FT of the two-point correlation function (Jing & Borner 2004), should be used to determine 
P(k) on small scales. 



4. SUMMARY 

In this paper, we have derived for the first time the formula [Eq. (17)] which relates the 
raw power spectrum (\5^(k)\ 2 ) estimated with FFT to the true power spectrum P(k). The 
formula shows clearly how the mass assignment modifies the power spectrum. The convo- 
lution of the density field with an assignment function W(r) reshapes the power spectrum 
(including the shot noise spectrum) by multiplying the factor |VF(k)| 2 ; the finite sampling 
of the convolved density field leads to the alias sum. We have described how to reconstruct 
P(k) from (\5 f (k)\ 2 ). For the NGP, CIC and TSC assignment functions, the shot noise 
-D 2 (k) can be expressed by simple analytical functions, therefore the shot noise can be easily 
corrected. To extract P(k) from the alias sum Y2 n W 2 (k + 2k^n.)P('k + 2/c^rn), we propose 
an iterative method. The method has been tested by applying it to an N-body simulation 
sample. Using different numbers of grid points, we have shown that the method can very 
successfully recover P(k) for all the NGP, CIC and TSC assignment functions. 
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Fig. 1. — The shot noise D 2 (k) estimated from 10 samples of Poisson distributed random 
points. Each symbol represents the result for each mass assignment, as indicated in the 
figure. The upper panel shows the result with the function {D 2 Ck)N/C^(]<.))d, i.e. the D 2 (k) 
scaled to l/iVCi(k). The estimated result agrees quite well with the analytical prediction 
(D 2 (\i)N/Cs{'k))d — 1- The lower panel compares the estimated (-D 2 (k)iV)d with our analyt- 
ical predictions for the NGP (solid line), CIC (dotted line) and TSC (dashed line) assignment 
functions. For CIC and TSC, we have used the approximate formulae of Eq. (21). 
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Fig. 2. — The six power spectra which we measured for an N-body simulation sample using 
three mass assignments (NGP, CIC and TSC) and two grids (64 3 and 256 3 grid points), k^f 
is the Nyquist wavenumber of 64 3 grid points. The upper panel shows the raw power spectra 
{\8f(k)\ 2 ) estimated directly from the FFT (Eq. 8). The lower panel shows the true power 
spectra P(k) reconstructed from the {\5- f (k)\ 2 } following the procedure described in the text. 
The six reconstructed P(k) agree so well that their curves overlay each other. 



